knitr::opts_knit$set(root.dir = '~/Documents/bgs_sims/')
knitr::opts_chunk$set(fig.width=16)
library(tidyverse)
library(ggridges)
library(viridis)
library(cowplot)
library(wesanderson)
orig_theme <- theme_set(theme_cowplot(font_size=20))
read_sim<-function(pidata,xidata,tajdata,roh=8.2e-10,neutral=FALSE,windows=20){
#read data
pi<-as.tibble(read.csv(pidata,header=T))
print(nrow(pi))
pi<-mutate(pi,stat=rep("pi",length(pi[,1])))
xi<-as.tibble(read.csv(xidata,header=T))
xi<-mutate(xi,stat=rep("xi",length(xi[,1])))
tajD <- as.tibble(read.csv(tajdata,header=T))
tajD<-mutate(tajD,stat=rep("tajD",length(tajD[,1])))
#munge into single tibble of means
data<-rbind(xi,pi,tajD) %>%
`colnames<-`(c("gen", 1:windows, "sim", "stat")) %>%
gather(window,value,-gen,-stat,-sim) %>%
group_by(gen,stat,window) %>%
summarize(mean_val=mean(value,na.rm = T)) %>%
mutate(centimorgen=(as.numeric(window)-1)*10000*roh*100 ,mean_val=mean_val)
return(data)
}
#neutral=TRUE
#pidata <- "results/tennessen_neutral_pi.csv"
#xidata <- "results/tennessen_neutral_xi.csv"
#tajdata <- "results/tennessen_neutral_tajD.csv"
Merge neutral and selected values
merge_sel_neutral <- function(sel_sim_data,neut_sim_data){
merged <- merge(sel_sim_data,neut_sim_data,by=c('gen','stat','window','centimorgen'))
merged <- mutate(merged,mean_val = if_else(stat=='tajD',mean_val.x-mean_val.y,mean_val.x/mean_val.y))
names(merged) <- c('gen', 'stat', 'window', 'centimorgen','bgs','neutral','mean_val')
merged
}
Plot selection landscape
simplot<-function(somedata,mystat,label){
ggplot(filter(somedata,stat==mystat), aes(x=centimorgen, y=mean_val, group = gen,color=gen)) +
scale_color_viridis(name='Generation') +
background_grid(major = "xy", minor = "none") +
# geom_smooth(se=F) +
geom_line() +
ylab(label) + xlab("Distance (cM)") #+
# scale_x_continuous(breaks = seq(-1.5,1.5,0.5))
}
simplot_grid <- function(somedata,events,event_names,ratio=TRUE){
if (ratio == TRUE){
pi_lab <- expression(bar(pi)/bar(pi)[neut])
xi_lab <- expression(bar(xi)/bar(xi)[neut])
d_lab <- expression(bar(D)-bar(D)[neut])
}
else{
pi_lab <- expression(bar(pi))
xi_lab <- expression(bar(xi))
d_lab <- expression(bar(D))
}
legend <- get_legend(simplot(somedata,"pi",expression(bar(pi)/bar(pi)[neut]))+
theme(legend.key.height= unit(2, "cm"),
legend.title = element_text(size=16),
legend.text = element_text(size=11)) +
scale_color_viridis(name='Generation',breaks = events, labels=event_names)
)
plot_grid(
plot_grid(
simplot(somedata,"pi",pi_lab)+theme(legend.position='none'),
simplot(somedata,"xi",xi_lab)+theme(legend.position="none"),
simplot(somedata,"tajD",d_lab)+theme(legend.position="none"),
ncol = 1,
align = 'hv'),
legend, rel_widths = c(85,15))
}
Function to plot ratios over time and by genetic distance
plot_ratios <- function(data,events,windows=20){
window_list <- c(1,5,10,20)
comb <- data %>%
filter(window%in%window_list,stat%in%c('pi','xi')) %>%
mutate(centimorgen=abs(centimorgen)) %>%
group_by(centimorgen,gen,stat) %>%
summarise(mean_val=mean(mean_val))
start_val <- comb %>%
filter(gen<0) %>%
group_by(stat,centimorgen) %>%
summarise(mean_val=mean(mean_val))
ggplot(comb,aes(gen,mean_val,color=stat)) +
geom_smooth(method = 'loess',span=0.1)+
# geom_line(size=1.5) +
geom_hline(data=data.frame(yint=start_val$mean_val,centimorgen=start_val$centimorgen),aes(yintercept =yint))+
labs(x=expression(paste('Generation (x ',N[anc],')')),y=expression(x/x[neut]),color='Statistics'
#,caption ='Summary statistics at different distances (cM) from selected region.\nVertical lines mark demographic events. Horizontal lines denote value before demographic changes'
) +
# xlim(-0.1,1.25) +
facet_grid(.~centimorgen) +
geom_vline(xintercept = events, size=0.5, alpha=0.2) +
scale_color_manual(values=c("#999999", "#E69F00"),
breaks=c("pi", "xi"),
labels=c(expression(bar(pi)),expression(bar(xi)))) +
theme(strip.background = element_rect(fill=alpha('lightblue',.2)),
strip.text.y = element_text(angle = 0,hjust = 0),
plot.caption=element_text(size=12,hjust = 0)) +
# scale_x_continuous(breaks=c(0,0.3,0.6,0.9,1.2)) +
NULL
}
Function to plot absolute values over time and by genetic distance
plot_over_time <- function(data,events,stats,windows=20){
center <- 1
window_list <- c(2,10,20)
comb <- data %>%
filter(window%in%window_list,stat==stats) %>%
mutate(centimorgen=abs(centimorgen)) %>%
group_by(centimorgen,gen,stat) %>%
summarise(neutral=mean(neutral),bgs=mean(bgs)) %>%
gather(selection,value,bgs,neutral)
start_val <- comb %>%
filter(gen<0) %>%
group_by(stat,centimorgen,selection) %>%
summarise(value=mean(value))
stats_label <- c(expression(bar(stats)))
ggplot(comb,aes(gen,value,color=selection)) +
geom_line(size=1.5) +
geom_hline(data=data.frame(yint=start_val$value,centimorgen=start_val$centimorgen,sel=start_val$selection),aes(yintercept =yint,color=sel))+
labs(x=expression(paste('Generation (x ',N[anc],')')),y=stats,color='Selection')+
facet_grid(.~centimorgen) +
geom_vline(xintercept = events, size=0.5, alpha=0.2) +
scale_color_manual(values=alpha(c("#ef8a62", "#67a9cf"),.9),
breaks=c("neutral", "bgs"),
labels=c('Neutral','BGS')) +
scale_x_continuous(breaks = seq(0,2,0.5))+
theme(strip.background = element_rect(fill=alpha('lightblue',.2)),
strip.text.y = element_text(angle = 0,hjust = 0),
plot.caption=element_text(size=12,hjust = 0))
}
demog_torres <- read.csv('demographies/torres.csv')
demog_torres <- mutate(demog_torres,model='torres')
demog_ten <- read.csv('demographies/tennessen.csv')
demog_ten <- mutate(demog_ten,model='tennessen')
demog_maize <- read.csv('demographies/maize_018.csv')
demog_maize <- mutate(demog_maize,model='maize')
demog <- rbind(demog_torres,demog_ten,demog_maize)
# Torres events
anc_ext <- 0
OOA_torres <- 0.437017*2
euro_bneck_torres <- OOA_torres + 0.109481*2
torres_events <- c(anc_ext,OOA_torres,euro_bneck_torres)
torres_event_names <- c('Ancestral expansion','OoA bottleneck','Euro bottleneck')
torres_event_sizes <- c(2.10709*18449,2.10709*18449,0.322295*18449)
torres_event_end<- c(1e5,1e5,2e3)
# Tennessen events
OOA_ten <- 0.264196*2
euro_bneck_ten <-OOA_ten+ 0.0763702*2
finel_ext <-euro_bneck_ten+ 0.0502841*2
tennessen_events <- c(anc_ext,OOA_ten,euro_bneck_ten,finel_ext)
tennessen_event_sizes <- c(1.98*7310,1.98*7310,0.254609*7310,1.339863*7310)
tennessen_event_end<- c(1e5,1e5,1e5,1e3)
tennessen_event_names <- c('Ancestral expansion','OoA bottleneck','Euro bottleneck','Final growth')
# maize events
maize_events <- c(0,0.12)
maize_event_end<- c(12278*0.04,12278*0.5)
maize_event_sizes <- c(12278,12278*2.98)
maize_event_names <- c('Domestication bottleneck','Maize present')
x <- filter(demog_maize,generation>0)
fit <- lm(log(x$N)~(x$generation))
x$pred <- exp(fit$coefficients[1]+x$generation * fit$coefficients[2])
ggplot()+
geom_line(data=x,aes(x=generation,y=pred),size=2) +
geom_line(data=x,aes(x=generation,y=N),color='red')
exp(fit$coefficients[1]+0.18 * fit$coefficients[2])
ggplot(demog,aes(x=generation,y=N,color=model))+geom_line(size=2) +
scale_y_log10()+
scale_color_discrete(breaks=c("tennessen", "torres","maize"),labels=c('Tennessen et al.','Torres et al.','Bessinger et al.'))+
# annotate tennessen
annotate("segment", x = tennessen_events, xend = tennessen_events,
y = tennessen_event_end, yend = tennessen_event_sizes,
size=1,alpha=0.6,arrow=arrow()) +
annotate("text", x = tennessen_events, y = tennessen_event_end,
label=tennessen_event_names,
angle=c(90,90,90,0),hjust=0,vjust=c(0.5,0.5,0.5,1)) +
# annotate torres
annotate("segment", x = torres_events, xend = torres_events,
y = torres_event_end, yend = torres_event_sizes,
size=1,alpha=0.6,arrow=arrow()) +
annotate("text", x = torres_events, y = torres_event_end,
label=torres_event_names,
angle=c(90,90,0),hjust=0,vjust=c(0.5,0.5,1)) +
# annotate maize
annotate("segment", x = maize_events, xend = maize_events,
y = maize_event_end, yend = maize_event_sizes,
size=1,alpha=0.6,arrow=arrow()) +
annotate("text", x = maize_events, y = maize_event_end,
label=maize_event_names,
angle=0,hjust=0,vjust=1) +
labs(x=expression(paste('Generation (x ',N[anc],')')),y= "Population size",color='Demography')
ggsave('figures/annotated_demographies.pdf',width = 16)
## Saving 16 x 5 in image
tennessen_neutral <-read_sim("results/tennessen_neutral_pi.csv",
"results/tennessen_neutral_xi.csv",
"results/tennessen_neutral_tajD.csv",
neutral=TRUE)
## [1] 610000
simplot_grid(tennessen_neutral,tennessen_events,tennessen_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
tennessen_bgs <-read_sim("results/tennessen_bgs_pi.csv",
"results/tennessen_bgs_xi.csv",
"results/tennessen_bgs_tajD.csv",
neutral=FALSE)
## [1] 610000
simplot_grid(tennessen_bgs,tennessen_events,tennessen_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
tennessen <- merge_sel_neutral(tennessen_bgs,tennessen_neutral)
legend <- get_legend(plot_over_time(tennessen,tennessen_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(tennessen,tennessen_events,'pi')+theme(legend.position='none'),
plot_over_time(tennessen,tennessen_events,'xi')+theme(legend.position='none'),
plot_over_time(tennessen,tennessen_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)
ggsave('figures/origVal_distance_summary_tennessen.pdf',width = 16,height = 12)
simplot_grid(tennessen,tennessen_events,tennessen_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
ggsave('figures/region_summary_tennessen.pdf',width = 12,heigh=9)
plot_ratios(tennessen,tennessen_events)
ggsave('figures/smooth_distance_summary_tennessen.pdf',width = 16,height = 7)
torres_neutral <-read_sim("results/torres_neutral_pi.csv",
"results/torres_neutral_xi.csv",
"results/torres_neutral_tajD.csv",
neutral=TRUE)
## [1] 2230000
simplot_grid(torres_neutral,torres_events,torres_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
torres_bgs <-read_sim("results/torres_bgs_pi.csv",
"results/torres_bgs_xi.csv",
"results/torres_bgs_tajD.csv",
neutral=FALSE)
## [1] 2230000
simplot_grid(torres_bgs,torres_events,torres_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
torres <- merge_sel_neutral(torres_bgs,torres_neutral)
legend <- get_legend(plot_over_time(torres,torres_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(torres,torres_events,'pi')+theme(legend.position='none'),
plot_over_time(torres,torres_events,'xi')+theme(legend.position='none'),
plot_over_time(torres,torres_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)
ggsave('figures/origVal_distance_summary_torres.pdf',width = 16,height = 12)
simplot_grid(torres,torres_events,torres_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
ggsave('figures/region_summary_torres.pdf',width = 12,heigh=9)
plot_ratios(torres,torres_events)
ggsave('figures/smooth_distance_summary_torres.pdf',width = 16,height=7)
maize_neutral <-read_sim("results/maize_neutral_pi.csv",
"results/maize_neutral_xi.csv",
"results/maize_neutral_tajD.csv",roh = 8.2e-9,
neutral=TRUE)
## [1] 165000
simplot_grid(maize_neutral,maize_events,maize_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
maize_bgs <-read_sim("results/maize_bgs_pi.csv",
"results/maize_bgs_xi.csv",
"results/maize_bgs_tajD.csv",roh = 8.2e-9,
neutral=FALSE)
## [1] 165000
simplot_grid(maize_bgs,maize_events,maize_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
maize <- merge_sel_neutral(maize_bgs,maize_neutral)
legend <- get_legend(plot_over_time(maize,maize_events,'pi',windows = 21))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(maize,maize_events,'pi',windows = 21)+theme(legend.position='none'),
plot_over_time(maize,maize_events,'xi',windows = 21)+theme(legend.position='none'),
plot_over_time(maize,maize_events,'tajD',windows = 21)+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)
ggsave('figures/origVal_distance_summary_maize.pdf',width = 16,height = 12)
simplot_grid(maize,maize_events,maize_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
ggsave('figures/region_summary_maize.pdf',width = 12,heigh=9)
plot_ratios(maize,maize_events,windows = 21)
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
ggsave('figures/smooth_distance_summary_maize.pdf',width = 16,height=7)
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
demog <- read.csv('demographies/generic_models.csv')
demog <- demog %>%
gather(model,size,Model.1,Model.2,Model.3,Model.4,Model.5,Model.6,Model.7,Model.8,Model.9,Model.10,Model.11,Model.12) %>%
mutate(bneck_length=case_when(model%in%c('Model.2','Model.3','Model.6','Model.7')~'1_Short',
model%in%c('Model.4','Model.5','Model.8','Model.9','Model.10','Model.11','Model.12')~'2_Long',
model=='Model.1'~'Constant'),
bneck_strength=case_when(model%in%c('Model.2','Model.4','Model.6','Model.8')~'1_Weak',
model%in%c('Model.3','Model.5','Model.7','Model.9')~'2_Strong',
model%in%c('Model.10','Model.11','Model.12')~'3_constant',
model=='Model.1'~'Constant'),
bneck_age = case_when(model%in%c('Model.6','Model.7','Model.8','Model.9')~'2_Late',
model%in%c('Model.2','Model.3','Model.4','Model.5','Model.10','Model.11','Model.12')~'1_Early',
model=='Model.1'~'Constant'),
size=size)
demog <- filter(demog,model!='Model.1')
length_labels <- c('1_Short' = 'Short','2_Long'='Long')
strength_labels <- c('1_Weak'= 'Weak','2_Strong'='Strong','3_constant'='Constant')
age_labels <- c('1_Early'='Old','2_Late'='Recent')
ggplot(demog,aes(generation,size,color=model))+
geom_hline(yintercept = 20000,color='grey80')+
geom_line(size=1.5) +
# geom_hline(yintercept = 400)+
scale_y_log10() +
theme(axis.text.x = element_text(angle = 45,hjust = 1,vjust = 1,face = 'bold')) +
labs(x=expression(paste('Generation (x ',N[anc],')')),
y='Population size',
color='Model')+
facet_grid(bneck_age~bneck_length+bneck_strength,space = 'free',
labeller = labeller(bneck_length=as_labeller(length_labels),
bneck_age = as_labeller(age_labels),
bneck_strength = as_labeller(strength_labels)))+
theme(strip.background = element_blank(), strip.placement = "outside")
## Warning: Removed 72000 rows containing missing values (geom_path).
ggsave('figures/generic_models/generic_demographies.pdf',width = 12)
## Saving 12 x 5 in image
## Warning: Removed 72000 rows containing missing values (geom_path).
model1_events <- c(0,0.9)
model1_event_names <- c(0,0.9)
model1_neutral <-read_sim("results/generic/model1_neutral_pi.csv",
"results/generic/model1_neutral_xi.csv",
"results/generic/model1_neutral_tajD.csv",
neutral=TRUE)
## [1] 2020000
simplot_grid(model1_neutral,model1_events,model1_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model1_bgs <-read_sim("results/generic/model1_bgs_pi.csv",
"results/generic/model1_bgs_xi.csv",
"results/generic/model1_bgs_tajD.csv",
neutral=FALSE)
## [1] 2020000
simplot_grid(model1_bgs,model1_events,model1_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model1 <- merge_sel_neutral(model1_bgs,model1_neutral)
head(model1)
## gen stat window centimorgen bgs neutral mean_val
## 1 -0.002 pi 1 0.00000 7.746545 13.31417 0.5818270
## 2 -0.002 pi 10 0.00738 8.398654 13.22437 0.6350891
## 3 -0.002 pi 11 0.00820 8.453766 13.27309 0.6369101
## 4 -0.002 pi 12 0.00902 8.593112 13.23504 0.6492697
## 5 -0.002 pi 13 0.00984 8.593202 13.17840 0.6520674
## 6 -0.002 pi 14 0.01066 8.640976 13.17171 0.6560254
var_pi <- function(theta,n){
my_var <- (((n+1)*theta)/(3*(n-1)))+((2*(n^2+n+3)*theta^2)/(9*n*(n-1)))
}
model1_pi_var <- var_pi(13.28,400)
sqrt(model1_pi_var)/sqrt(1000)
## [1] 0.2093724
theta <- 13.28
n=400
neutral_pi <- read.csv("results/generic/model1_neutral_pi.csv")
head(neutral_pi)
## gen X0 X1 X2 X3 X4 X5 X6 X7 X8 X9
## 1 -0.008 8.130 6.455 7.739 4.571 6.413 9.515 11.278 23.822 27.641 22.982
## 2 -0.005 7.494 6.346 7.435 4.393 5.506 8.245 10.167 23.088 27.386 22.314
## 3 -0.002 7.480 6.145 7.494 4.478 5.670 8.304 9.948 22.789 26.858 22.266
## 4 0.000 7.707 5.784 6.982 4.600 6.162 8.587 10.491 23.800 27.238 22.731
## 5 0.002 7.785 6.102 7.222 4.618 6.089 8.676 10.636 24.006 27.406 23.177
## 6 0.005 7.442 6.119 7.363 4.299 5.653 8.100 10.081 23.655 27.590 23.674
## X10 X11 X12 X13 X14 X15 X16 X17 X18 X19
## 1 18.017 15.643 17.159 19.201 32.965 24.574 18.818 12.661 10.248 12.846
## 2 18.403 15.911 18.317 19.211 33.172 24.907 20.046 13.443 10.486 13.373
## 3 18.228 16.044 17.707 19.227 32.763 24.874 20.149 13.378 10.799 13.164
## 4 18.238 15.721 18.088 19.238 32.783 24.502 19.329 13.527 10.915 13.006
## 5 18.014 15.779 17.278 19.180 32.910 24.607 18.657 13.024 10.842 12.909
## 6 18.195 15.877 17.559 18.757 32.426 24.228 19.491 13.388 10.862 13.215
## replicate
## 1 2511
## 2 2511
## 3 2511
## 4 2511
## 5 2511
## 6 2511
tail(neutral_pi)
## gen X0 X1 X2 X3 X4 X5 X6 X7
## 2019995 0.988 17.257 11.415 11.729 10.725 10.443 8.615 13.411 17.169
## 2019996 0.990 17.136 11.437 11.878 10.838 10.095 7.859 12.341 15.283
## 2019997 0.992 17.798 11.555 11.609 10.771 10.732 9.078 14.224 17.401
## 2019998 0.995 18.433 12.233 12.231 11.150 10.783 8.834 14.301 17.977
## 2019999 0.998 18.240 11.904 11.882 10.980 11.125 8.966 14.447 18.557
## 2020000 1.000 18.560 12.118 12.088 11.218 11.242 9.199 14.752 19.128
## X8 X9 X10 X11 X12 X13 X14 X15 X16
## 2019995 18.808 16.767 15.301 11.958 16.248 12.965 11.644 10.300 12.018
## 2019996 17.854 15.532 14.030 11.308 14.759 11.708 10.579 9.313 10.955
## 2019997 19.322 17.499 16.009 12.349 15.984 13.000 11.906 10.267 11.816
## 2019998 19.395 17.315 15.766 12.227 16.132 12.092 10.855 9.599 11.235
## 2019999 19.688 17.851 16.180 12.410 16.128 12.221 11.489 9.981 11.905
## 2020000 19.856 17.729 16.685 12.892 17.316 13.248 12.467 10.995 12.573
## X17 X18 X19 replicate
## 2019995 11.199 8.521 6.381 4308
## 2019996 11.025 8.197 6.068 4308
## 2019997 12.624 9.186 6.088 4308
## 2019998 11.467 8.690 6.455 4308
## 2019999 11.533 8.802 6.414 4308
## 2020000 11.824 9.311 6.627 4308
neutral_pi %>%
gather(value,window,-replicate,-gen) %>%
group_by(replicate,value) %>%
summarise(window=mean(window))%>%
filter(window>13.28+sqrt(model1_pi_var)*2)
## # A tibble: 3,060 x 3
## # Groups: replicate [1,460]
## replicate value window
## <int> <chr> <dbl>
## 1 2 X3 31.6
## 2 2 X6 26.8
## 3 4 X10 27.5
## 4 4 X9 29.4
## 5 10 X10 28.0
## 6 10 X4 33.0
## 7 10 X5 29.8
## 8 10 X6 27.6
## 9 10 X7 34.4
## 10 13 X15 33.7
## # ... with 3,050 more rows
3060/100000
## [1] 0.0306
length(unique(neutral_pi$replicate))
## [1] 5000
ggplot(neutral_pi,aes(gen,X1,group=replicate))+geom_line()+
geom_ribbon(aes(ymax = 13.28+sqrt(model1_pi_var)*2, ymin = 13.28-(sqrt(model1_pi_var)*2)))
neutral_pi %>%
group_by(replicate)%>%
summarise_at(c("X1", "X2"), mean, na.rm = TRUE) %>%
filter(X1 >13.28+sqrt(model1_pi_var)*2, X2 >13.28+sqrt(model1_pi_var)*2) %>%
count()
## # A tibble: 1 x 1
## n
## <int>
## 1 72
72/5000
## [1] 0.0144
mean(neutral_pi$window)
## Warning in mean.default(neutral_pi$window): argument is not numeric or
## logical: returning NA
## [1] NA
legend <- get_legend(plot_over_time(model1,model1_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model1,model1_events,'pi')+theme(legend.position='none')+geom_hline(yintercept = 13.28),
plot_over_time(model1,model1_events,'xi')+theme(legend.position='none')+geom_hline(yintercept = 13.28),
plot_over_time(model1,model1_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)
ggsave('figures/generic_models/origVal_distance_summary_model1.pdf',width = 16,height = 12)
ggplot(filter(model1,stat=='pi',window%in%c(1,10,20)),aes(gen,neutral,group=window))+
geom_ribbon(aes(ymax = 13.28+sqrt(model1_pi_var), ymin = 13.28-sqrt(model1_pi_var)))+
geom_ribbon(aes(ymax = neutral + sd_neutral/sqrt(1000), ymin = neutral - sd_neutral/sqrt(1000),fill=as.factor(as.numeric(window))),alpha=.2) +geom_line()+
geom_hline(yintercept = 13.28)+
labs(fill='window ')
ggsave('test/origVal_distance_summary_model1.pdf',width = 10,height = 10)
simplot_grid(model1,model1_events,model1_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
ggsave('figures/generic_models/region_summary_model1.pdf',width = 12,heigh=9)
plot_ratios(model1,model1_events)
ggsave('figures/generic_models/smooth_distance_summary_model1.pdf',width = 16,height=7)
model2_events <- c(0,0.5)
model2_event_names <- c('Bottleneck','Recovered')
model2_neutral <-read_sim("results/generic/model2_neutral_pi.csv",
"results/generic/model2_neutral_xi.csv",
"results/generic/model2_neutral_tajD.csv",
neutral=TRUE)
## [1] 2020000
simplot_grid(model2_neutral,model2_events,model2_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model2_bgs <-read_sim("results/generic/model2_bgs_pi.csv",
"results/generic/model2_bgs_xi.csv",
"results/generic/model2_bgs_tajD.csv",
neutral=FALSE)
## [1] 2020000
simplot_grid(model2_bgs,model2_events,model2_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model2 <- merge_sel_neutral(model2_bgs,model2_neutral)
legend <- get_legend(plot_over_time(model2,model2_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model2,model2_events,'pi')+theme(legend.position='none'),
plot_over_time(model2,model2_events,'xi')+theme(legend.position='none'),
plot_over_time(model2,model2_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)
ggsave('figures/generic_models/origVal_distance_summary_model2.pdf',width = 16,height = 12)
simplot_grid(model2,model2_events,model2_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
ggsave('figures/generic_models/region_summary_model2.pdf',width = 12,heigh=9)
plot_ratios(model2,model2_events)
ggsave('figures/generic_models/smooth_distance_summary_model2.pdf',width = 16,height=7)
model3_events <- c(0,0.6667)
model3_event_names <- c('Bottleneck','Recovered')
model3_neutral <-read_sim("results/generic/model3_neutral_pi.csv",
"results/generic/model3_neutral_xi.csv",
"results/generic/model3_neutral_tajD.csv",
neutral=TRUE)
## [1] 2020000
simplot_grid(model3_neutral,model3_events,model3_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model3_bgs <-read_sim("results/generic/model3_bgs_pi.csv",
"results/generic/model3_bgs_xi.csv",
"results/generic/model3_bgs_tajD.csv",
neutral=FALSE)
## [1] 2020000
simplot_grid(model3_bgs,model3_events,model3_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model3 <- merge_sel_neutral(model3_bgs,model3_neutral)
legend <- get_legend(plot_over_time(model3,model3_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model3,model3_events,'pi')+theme(legend.position='none'),
plot_over_time(model3,model3_events,'xi')+theme(legend.position='none'),
plot_over_time(model3,model3_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)
ggsave('figures/generic_models/origVal_distance_summary_model3.pdf',width = 16,height = 12)
simplot_grid(model3,model3_events,model3_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
ggsave('figures/generic_models/region_summary_model3.pdf',width = 12,heigh=9)
plot_ratios(model3,model3_events)
ggsave('figures/generic_models/smooth_distance_summary_model3.pdf',width = 16,height=7)
model4_events <- c(0,0.05,0.5250)
model4_event_names <- c('Bottleneck','Growth','Recovered')
model4_neutral <-read_sim("results/generic/model4_neutral_pi.csv",
"results/generic/model4_neutral_xi.csv",
"results/generic/model4_neutral_tajD.csv",
neutral=TRUE)
## [1] 2020000
simplot_grid(model4_neutral,model4_events,model4_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model4_bgs <-read_sim("results/generic/model4_bgs_pi.csv",
"results/generic/model4_bgs_xi.csv",
"results/generic/model4_bgs_tajD.csv",
neutral=FALSE)
## [1] 2020000
simplot_grid(model4_bgs,model4_events,model4_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model4 <- merge_sel_neutral(model4_bgs,model4_neutral)
legend <- get_legend(plot_over_time(model4,model4_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model4,model4_events,'pi')+theme(legend.position='none'),
plot_over_time(model4,model4_events,'xi')+theme(legend.position='none'),
plot_over_time(model4,model4_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)
ggsave('figures/generic_models/origVal_distance_summary_model4.pdf',width = 16,height = 12)
simplot_grid(model4,model4_events,model4_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
ggsave('figures/generic_models/region_summary_model4.pdf',width = 12,heigh=9)
plot_ratios(model4,model4_events)
ggsave('figures/generic_models/smooth_distance_summary_model4.pdf',width = 16,height=7)
model5_events <- c(0,0.05,0.6833)
model5_event_names <- c('Bottleneck','Growth','Recovered')
model5_neutral <-read_sim("results/generic/model5_neutral_pi.csv",
"results/generic/model5_neutral_xi.csv",
"results/generic/model5_neutral_tajD.csv",
neutral=TRUE)
## [1] 2020000
simplot_grid(model5_neutral,model5_events,model5_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model5_bgs <-read_sim("results/generic/model5_bgs_pi.csv",
"results/generic/model5_bgs_xi.csv",
"results/generic/model5_bgs_tajD.csv",
neutral=FALSE)
## [1] 2020000
simplot_grid(model5_bgs,model5_events,model5_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model5 <- merge_sel_neutral(model5_bgs,model5_neutral)
legend <- get_legend(plot_over_time(model5,model5_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model5,model5_events,'pi')+theme(legend.position='none'),
plot_over_time(model5,model5_events,'xi')+theme(legend.position='none'),
plot_over_time(model5,model5_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)
ggsave('figures/generic_models/origVal_distance_summary_model5.pdf',width = 16,height = 12)
simplot_grid(model5,model5_events,model5_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
ggsave('figures/generic_models/region_summary_model5.pdf',width = 12,heigh=9)
plot_ratios(model5,model5_events)
ggsave('figures/generic_models/smooth_distance_summary_model5.pdf',width = 16,height=7)
model6_events <- c(0.,0.05)
model6_event_names <- c('Bottleneck','Recovered')
model6_neutral <-read_sim("results/generic/model6_neutral_pi.csv",
"results/generic/model6_neutral_xi.csv",
"results/generic/model6_neutral_tajD.csv",
neutral=TRUE)
## [1] 220000
simplot_grid(model6_neutral,model6_events,model6_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model6_bgs <-read_sim("results/generic/model6_bgs_pi.csv",
"results/generic/model6_bgs_xi.csv",
"results/generic/model6_bgs_tajD.csv",
neutral=FALSE)
## [1] 220000
simplot_grid(model6_bgs,model6_events,model6_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model6 <- merge_sel_neutral(model6_bgs,model6_neutral)
legend <- get_legend(plot_over_time(model6,model6_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model6,model6_events,'pi')+theme(legend.position='none'),
plot_over_time(model6,model6_events,'xi')+theme(legend.position='none'),
plot_over_time(model6,model6_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)
ggsave('figures/generic_models/origVal_distance_summary_model6.pdf',width = 16,height = 12)
simplot_grid(model6,model6_events,model6_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
ggsave('figures/generic_models/region_summary_model6.pdf',width = 12,heigh=9)
plot_ratios(model6,model6_events)
ggsave('figures/generic_models/smooth_distance_summary_model6.pdf',width = 16,height=7)
model7_events <- c(0.0,0.0667)
model7_event_names <- c('Bottleneck','Recovered')
model7_neutral <-read_sim("results/generic/model7_neutral_pi.csv",
"results/generic/model7_neutral_xi.csv",
"results/generic/model7_neutral_tajD.csv",
neutral=TRUE)
## [1] 220000
simplot_grid(model7_neutral,model7_events,model7_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model7_bgs <-read_sim("results/generic/model7_bgs_pi.csv",
"results/generic/model7_bgs_xi.csv",
"results/generic/model7_bgs_tajD.csv",
neutral=FALSE)
## [1] 220000
simplot_grid(model7_bgs,model7_events,model7_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model7 <- merge_sel_neutral(model7_bgs,model7_neutral)
legend <- get_legend(plot_over_time(model7,model7_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model7,model7_events,'pi')+theme(legend.position='none'),
plot_over_time(model7,model7_events,'xi')+theme(legend.position='none'),
plot_over_time(model7,model7_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)
ggsave('figures/generic_models/origVal_distance_summary_model7.pdf',width = 16,height = 12)
simplot_grid(model7,model7_events,model7_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
ggsave('figures/generic_models/region_summary_model7.pdf',width = 12,heigh=9)
plot_ratios(model7,model7_events)
ggsave('figures/generic_models/smooth_distance_summary_model7.pdf',width = 16,height=7)
model8_events <- c(0.0,0.05,0.0750)
model8_event_names <- c('Bottleneck','Growth','Recovered')
model8_neutral <-read_sim("results/generic/model8_neutral_pi.csv",
"results/generic/model8_neutral_xi.csv",
"results/generic/model8_neutral_tajD.csv",
neutral=TRUE)
## [1] 220000
simplot_grid(model8_neutral,model8_events,model8_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model8_bgs <-read_sim("results/generic/model8_bgs_pi.csv",
"results/generic/model8_bgs_xi.csv",
"results/generic/model8_bgs_tajD.csv",
neutral=FALSE)
## [1] 220000
simplot_grid(model8_bgs,model8_events,model8_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model8 <- merge_sel_neutral(model8_bgs,model8_neutral)
legend <- get_legend(plot_over_time(model8,model8_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model8,model8_events,'pi')+theme(legend.position='none'),
plot_over_time(model8,model8_events,'xi')+theme(legend.position='none'),
plot_over_time(model8,model8_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)
ggsave('figures/generic_models/origVal_distance_summary_model8.pdf',width = 16,height = 12)
simplot_grid(model8,model8_events,model8_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
ggsave('figures/generic_models/region_summary_model8.pdf',width = 12,heigh=9)
plot_ratios(model8,model8_events)
ggsave('figures/generic_models/smooth_distance_summary_model8.pdf',width = 16,height=7)
model9_events <- c(0.0,0.05,0.0833)
model9_event_names <- c('Bottleneck','Growth','Recovered')
model9_neutral <-read_sim("results/generic/model9_neutral_pi.csv",
"results/generic/model9_neutral_xi.csv",
"results/generic/model9_neutral_tajD.csv",
neutral=TRUE)
## [1] 220000
simplot_grid(model9_neutral,model9_events,model9_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model9_bgs <-read_sim("results/generic/model9_bgs_pi.csv",
"results/generic/model9_bgs_xi.csv",
"results/generic/model9_bgs_tajD.csv",
neutral=FALSE)
## [1] 220000
simplot_grid(model9_bgs,model9_events,model9_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model9 <- merge_sel_neutral(model9_bgs,model9_neutral)
legend <- get_legend(plot_over_time(model9,model9_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model9,model9_events,'pi')+theme(legend.position='none'),
plot_over_time(model9,model9_events,'xi')+theme(legend.position='none'),
plot_over_time(model9,model9_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)
ggsave('figures/generic_models/origVal_distance_summary_model9.pdf',width = 16,height = 12)
simplot_grid(model9,model9_events,model9_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
ggsave('figures/generic_models/region_summary_model9.pdf',width = 12,heigh=9)
plot_ratios(model9,model9_events)
ggsave('figures/generic_models/smooth_distance_summary_model9.pdf',width = 16,height=7)
model10_events <- c(0.0)
model10_event_names <- c('Bottleneck')
model10_neutral <-read_sim("results/generic/model10_neutral_pi.csv",
"results/generic/model10_neutral_xi.csv",
"results/generic/model10_neutral_tajD.csv",
neutral=TRUE)
## [1] 2020000
simplot_grid(model10_neutral,model10_events,model10_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model10_bgs <-read_sim("results/generic/model10_bgs_pi.csv",
"results/generic/model10_bgs_xi.csv",
"results/generic/model10_bgs_tajD.csv",
neutral=FALSE)
## [1] 2020000
simplot_grid(model10_bgs,model10_events,model10_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model10 <- merge_sel_neutral(model10_bgs,model10_neutral)
legend <- get_legend(plot_over_time(model10,model10_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model10,model10_events,'pi')+theme(legend.position='none'),
plot_over_time(model10,model10_events,'xi')+theme(legend.position='none'),
plot_over_time(model10,model10_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)
ggsave('figures/generic_models/origVal_distance_summary_model10.pdf',width = 16,height = 12)
simplot_grid(model10,model10_events,model10_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
ggsave('figures/generic_models/region_summary_model10.pdf',width = 12,heigh=9)
plot_ratios(model10,model10_events)
ggsave('figures/generic_models/smooth_distance_summary_model10.pdf',width = 16,height=7)
model11_events <- c(0.0)
model11_event_names <- c('Bottleneck')
model11_neutral <-read_sim("results/generic/model11_neutral_pi.csv",
"results/generic/model11_neutral_xi.csv",
"results/generic/model11_neutral_tajD.csv",
neutral=TRUE)
## [1] 2020000
simplot_grid(model11_neutral,model11_events,model11_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model11_bgs <-read_sim("results/generic/model11_bgs_pi.csv",
"results/generic/model11_bgs_xi.csv",
"results/generic/model11_bgs_tajD.csv",
neutral=FALSE)
## [1] 2020000
simplot_grid(model11_bgs,model11_events,model11_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model11 <- merge_sel_neutral(model11_bgs,model11_neutral)
legend <- get_legend(plot_over_time(model11,model11_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model11,model11_events,'pi')+theme(legend.position='none'),
plot_over_time(model11,model11_events,'xi')+theme(legend.position='none'),
plot_over_time(model11,model11_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)
ggsave('figures/generic_models/origVal_distance_summary_model11.pdf',width = 16,height = 12)
simplot_grid(model11,model11_events,model11_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
ggsave('figures/generic_models/region_summary_model11.pdf',width = 12,heigh=9)
plot_ratios(model11,model11_events)
ggsave('figures/generic_models/smooth_distance_summary_model11.pdf',width = 16,height=7)
model12_events <- c(0.0)
model12_event_names <- c('Growth')
model12_neutral <-read_sim("results/generic/model12_neutral_pi.csv",
"results/generic/model12_neutral_xi.csv",
"results/generic/model12_neutral_tajD.csv",
neutral=TRUE)
## [1] 2020000
simplot_grid(model12_neutral,model12_events,model12_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model12_bgs <-read_sim("results/generic/model12_bgs_pi.csv",
"results/generic/model12_bgs_xi.csv",
"results/generic/model12_bgs_tajD.csv",
neutral=FALSE)
## [1] 2020000
simplot_grid(model12_bgs,model12_events,model12_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
model12 <- merge_sel_neutral(model12_bgs,model12_neutral)
legend <- get_legend(plot_over_time(model12,model12_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model12,model12_events,'pi')+theme(legend.position='none'),
plot_over_time(model12,model12_events,'xi')+theme(legend.position='none'),
plot_over_time(model12,model12_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)
ggsave('figures/generic_models/origVal_distance_summary_model12.pdf',width = 16,height = 12)
simplot_grid(model12,model12_events,model12_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
ggsave('figures/generic_models/region_summary_model12.pdf',width = 12,heigh=9)
plot_ratios(model12,model12_events)
ggsave('figures/generic_models/smooth_distance_summary_model12.pdf',width = 16,height=7)